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A scalable computational approach to the simulation of propellant tank sloshing 
dynamics in microgravity is presented. In this work, we use the lattice Boltzmann 
equation (LBE) to approximate the behavior of two-phase, single-component isother- 
mal flows at very low Bond numbers. Through the use of a non-ideal gas equa- 
tion of state and a modified multiple relaxation time (MRT) collision operator, the 
proposed method can simulate thermodynamically consistent phase transitions at 
temperatures and density ratios consistent with typical spacecraft cryogenic pro- 
pellants, for example, liquid oxygen. Determination of the tank forces and mo- 
ments relies upon the global momentum conservation of the fluid domain, and 
a parametric wall wetting model allows tuning of the free surface contact angle. 
Development of the interface is implicit and no interface tracking approach is re- 
quired. Numerical examples illustrate the method’s application to predicting bulk 
fluid motion including lateral propellant slosh in low-g conditions. 


1 INTRODUCTION 

The modeling and predietion of the behavior of fluids in mierogravity eontinues to be a ehallenge 
in the design of spaeeeraft systems.^ In the mierogravity environment, hydrodynamie regimes can 
be described by the non-dimensional parameter of Bond number (Bo), characterizing the relative 
magnitude of the gravitational acceleration versus the capillary forces present in the liquid.^ At 
very low Bond numbers {Bo < 10), the hydrodynamics are dominated by the surface tension and a 
qualitative change in behavior is observed. Liquid free surface interfaces become characteristically 
curved, and most propellants approach a near zero contact angle with solid objects (such as tank 
walls). The dominant time scale of the liquid dynamics increases into the tens or hundreds of 
seconds, and characteristic flow velocities and the Mach number are very small (M <C 1). In these 
conditions, computational fluid mechanics (CFM) approaches are required to predict the motion of 
the bulk fluid mass and its effect on the spacecraft when displaced from equilibrium. 

The lattice Boltzmann equation (LBE) has recently emerged as a promising alternative to tradi- 
tional approaches to computational fluid mechanics. The LBE is a method by which the continuum 
fluid transport phenomena, i.e. the Navier-Stokes equations, can be approximated as a solution of a 
discretized nonlinear difference equation based upon the kinetic theory of gases. While the EBE is 
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typically restricted to very low- velocity flows, it does provide several unique advantages over tradi- 
tional solvers. First, the meshing of complex geometry is performed on a regular cartesian lattice of 
fluid cells, each having uniform volume in the fluid domain. As such, computations involving flux 
across the boundary of adjacent cells are considerably simplified. Secondly, LBE has the advantage 
of data locality; LBE-based flow solvers are not required to solve a global continuity equation at 
each time step. Einally, EBE is relatively simple to implement and computationally efficient. 

In the following sections, the development of an EBE-based flow solver for spacecraft propellant 
dynamics will be detailed. In Section 2, the theory and basic implementation details of the EBE 
will be introduced. In Section 3, the method of introducing multi-phase behavior into the EBE will 
be discussed. In Section 4, the results of test cases that compare the outputs of the EBE flow solver 
with theoretical predictions will be presented. In Section 5, a summary of the present research will 
be provided along with some opportunities for forward work. 

2 THE LATTICE BOLTZMANN METHOD 

The lattice Boltzmann equation (LBE) is a discretization of the continuous Boltzmann equation, 
describing particle dynamics on a molecular scale. The Boltzmann equation is given by 

||-+C^Vx/ + a^V^/ = F!(/) (1) 

where /(x, t) is the molecular velocity distribution function (DE) in the phase space (x, where 
X G is the spatial position and ^ G is the velocity. The derivation of the Boltzmann 
equation follows from the statistical kinetic theory of dilute gases in a closed domain. Here, a is the 
acceleration due to the applied body force at the spatial location x, and the collision integral is given 
by n. The right-hand side of the Boltzmann equation is the collision term describing the short-range 
molecular interactions of the velocity distributions assuming the modeled fluid is a dilute gas. 

The direct computation of the collision operator is, in general, intractable for the continuous 
Boltzmann equation. However, it can be approximated by a relaxation operation that preserves the 
hydrodynamic moments that are invariants of the collision. In the simplest models, such as the BGK 
(Bhatnagar-Gross-Krook) approximation, the collision is approximated by a linear relaxation to the 
equilibrium distribution, which is related to the temperature and velocity of the flow. 

The lattice Boltzmann equation follows from discretization of the 2Z)-dimensional phase space 
and the local approximation of the resultant linear system of ordinary differential equations in dis- 
crete time, and the approximation of the equilibrium distribution function consistent with that ve- 
locity discretization. Eor the present model D = 2 and a velocity discretization of 9 directions in 
two dimensions is chosen. The spatial discretization is applied on a regular lattice of size 6x with a 
temporal discretization 5t. The lattice structure shown in Eigure 1 is known as the D2Q9 model. 

The lattice parameter c = 5x/6t defines the characteristic speed associated with the velocity dis- 
cretization Oj, i = 0, 1, . . . 8. In this discretization, the cell at spatial location x^ is described 
by a distribution function f(x^, t) G M®, the velocity distribution function at the lattice site x^. The 
component of the discretized distribution function describes the density of particles at x^ having 
velocity e*. The approximation of the equilibrium distribution f®'! (namely, the Maxwell equilib- 
rium) is carried out such that the kinetic hydrodynamic moments are consistently approximated after 
discretization. The continuous Maxwell distribution is expanded to third order in the velocity; the 
truncated equilibrium approximation then converted into a discretized form using a Gauss-Hermite 
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Figure 1. D2Q9 Lattice Structure 


quadrature. The solution of the discretized Boltzmann equation can then be approximated by 

f(xfc + eJf,fo + <5f) = f(xfc,fo) - A(f(xfc,to) - f®'*(xfc,Ufc,p(xfc)))<5L (2) 

where is the local velocity and the collision matrix A has been introduced. The multiple scalar 
advection operations are implied by the operator e for notational convenience. The conserved hy- 
drodynamic moments are given by 


p(xfc) 

2 

(3) 

p(xfc)ufc 


(4) 


i 


and the speed of sound is Cg = cjy/3. 

Equation 2 is the lattice Boltzmann equation. Under certain conditions, it is possible to recover 
the macroscopic transport and continuity equations using the Chapman-Enskog expansion of the 
EBE.^ The approximation of incompressible flow is sensitive to velocity, partly due to the veloc- 
ity truncation in the series expansion of the equilibrium distribution function. This implies that 
compressibility error is the dominant error source in the application of EBE. 

The EBE is typically implemented in two steps comprising advection (streaming) and collision 
(relaxation) using two copies of the domain, f and f'. The boundary conditions are applied at lattice 
sites where fluid cells are adjacent to a prescribed boundary, such as a wall. In almost all numerical 
implementations, the units associated with the lattice are chosen such that 5x = 5t = c = 1 and 
the mean density p is on the order of 1 . This greatly reduces the computational burden and provides 
good numerical conditioning for the requisite computations. 

In the following developments, the spatial dependency of the hydrodynamic quantities will be 
denoted with a subscript k for brevity, and the temporal dependency is implied. 

2.1 Relaxation Operation 

In implementations of the EBE using the BGK approximation, the collision matrix is replaced by 
a scalar relaxation frequency ^ where r is the relaxation time. In this single relaxation time 
(SRT) model, all populations relax toward equilibrium at the same rate. It has been recognized that 


3 


while the SRT model is simple to implement, the relaxation of all populations at the same rate does 
not allow for preeise eontrol of the higher-order kinetie moments, whieh leads to poor numerieal 
stability and impreeise boundary eonditions.^’^ In order to overeome this limitation, the relaxation 
operation ean be performed in a moment spaee using a suitable linear transformation M. Stability 
ean be markedly improved by relaxing the moments at different rates. A set of relaxation rates Si 
are ehosen to satisfy eonservation eonstraints and maintain the neeessary symmetry.^ 

The hydrodynamie and kinetie moments ean be assembled in a veetor form and a linear trans- 
formation ean be derived sueh that m = Mf where M is the moment matrix. The ^th 

eollision 

operation is then performed in a moment spaee; 

A (ffc[n] - r'^(xfc,Ufc,/9fc)) = M“^S (rrifc[n] - Ufc, pfc)) . (5) 

Sinee M is a similarity transformation ehosen to orthogonalize the modes, S is a diagonal matrix 
whose eigenvalues s* are the relaxation parameters. In the isothermal LBE, the eonserved hydro- 
dynamie moments are only the density p and the momentum jx, jy\ thus the ehoiee of parameters 
So, S 3 , and S 5 are uneonstrained. The parameters S 4 = se, S 7 = sg must be equated due to symme- 
try and are related to the energy flux and kinematie viseosity, respeetively. Finally, there is a free 
ehoiee of parameters si and S 2 , whieh are related to the bulk viseosity and the energy. This seheme 
is known as the multiple relaxation time (MRT) model. 

Usually, the ehoiee of these parameters depends on a stability analysis eondueted through lin- 
earization of the eollision operator in order to determine the s* that provide the most stable eigen- 
values of the linearized LBE.^ Alternatively, they ean be seleeted to provide the best boundary 
eondition aeeuraey. It has been shown that a proper ehoiee of these parameters ean drastieally im- 
prove numerieal stability of the EBE while maintaining good physieal aeeuraey of the solution, 
allowing the simulation of very high Reynolds numbers in the single-phase ease.^’^ 

2.2 Body Forces 

The proper ineorporation of external forees into the EBE is an area of aetive research.’ Since the 
derivation of the EBE does not explicitly account for the incorporation of the acceleration term in 
Equation 1, various methods of incorporating these effects have been proposed. 

Kupershtokh* devised a method to incorporate the force that simultaneously acts on the advection 
step and the equilibrium distribution to ensure that the entire lattice, if in local equilibrium, remains 
in local equilibrium after a force is applied. In this Exact Difference Method (EDM), the shift in 
velocity due to the applied force is first computed (with 5t = 1) as Au^ = F^/p^, noting that 
Aufc = g everywhere if a uniform acceleration is present. The lattice distribution functions are first 
shifted by an amount equivalent to the change in equilibrium due to this velocity increment without 
advection. The velocity field is then shifted and the collision is performed using the modified 
velocity u^. 

The computational implication of the increase in accuracy is that the equilibrium distribution 
must be computed at least twice. However, these steps can be combined into a single collision 
equation; 


ffc[n + l] = (I-M ^SM) - ^‘^(xfc, Ufc, p^)) -h f®'^(xfc, u^, p^). 


( 6 ) 



In our model, this mechanization provides about a 20% performance improvement over the orig- 
inal three-step process. In particular, the matrix I— M“^SM may be precomputed, saving consid- 
erable computational expense at each lattice site. 

In an advanced MRT model, especially when used for multiphase flow, it is necessary to relax not 
only the hydrodynamic modes at different rates, but individual cells must have different effective 
relaxation frequencies Uk- If this is the case, the matrix A can be expanded as 

Aq - h cufcAAfc = (So + ASfc) M (7) 

where So is the bulk hydrodynamic part of the collision matrix, with S 7 = ss = 0. The matrix A A^ 
therefore relaxes only the diagonal and off-diagonal components of the stress tensor p^:x and p^y, 
respectively. Letting Affc[n] = f(.[n] — u^., p^), the MRT collision equation has the final 

form 


ik[n + l] = (I-Ao) Affc[n] + r'»(xfc,u^,p(xfc)) - cufcAAfcAffc[n]. (8) 

Global conservation of momentum can be used to deduce the net force on the domain, which is 
of primary importance for spacecraft dynamics applications. The lattice site momentum density can 
be used to compute the net momentum at time step n as p[n] = Ylk Pk’'^k\ since the fluid mass and 
time step are equal to unity in the LBE, the net change in momentum density is equal to the internal 
force from Newton’s law. The total domain force is therefore the difference of the internal force 
and the body force; 

k 


3 MULTI-PHASE FLOW 

The incorporation of multi-phase phenomena into the LBE is an emerging topic and has received 
much recent attention in the literature. Numerous methods for incorporating phase transitions in 
the EBE have been proposed.^ These can generally be divided into two major categories of explicit 
and implicit interface methods. In the explicit interface methods, similar to the approaches used 
for Navier-Stokes codes, the two fluid species are tracked via markers or coloring methods and 
the interface is defined at each time step. Either mass transfer across the interface is explicitly 
computed, or the interface tracking is performed with level sets. 

Implicit interface methods such as the Shan-Chen model^° simulate multiple phases uniformly 
in the lattice via the use of a pseudopotential function which models intramolecular interactions. 
A real gas equation of state can be used to elicit thermodynamically consistent phase transitions, 
which allows direct parameterization of the temperature.^^’ Effects such as phase separation, 
evaporation, and condensation occur continuously in the lattice while mass is conserved globally. 
This capability is particularly attractive for the modeling of cryogens, where the fluid temperature 
T <Tc may not be extremely subcritical and phase transitions must be accurately captured. 

The pseudopotential model is a diffuse interface model; that is, the interface appears in the lattice 
as a density gradient over a few cells or tens of cells. This increases the required minimum lattice 
resolution to obtain sharp interface definition if required by the application. A disadvantage of the 
pseudopotential model is the interdependency of the interface thickness and the surface tension pa- 
rameter. In addition, high pseudopotential forces are required to maintain the interface, degrading 
the numerical stability. Of course, the basic EBE is isothermal; convective effects, which can be 
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important in cryogenic flows in zero-g, are not captured.^ However, the simplicity of implementa- 
tion makes the pseudopotential model attractive in the present application insofar as it can be used 
to characterize bulk fluid behavior in response to transient accelerations. 


3.1 Pseudopotential Model 

In the formulation of a pseudopotential multiphase model based on the Shan-Chen (SC)**^ scheme, 
the intramolecular interaction force is related to a scalar function that depends on the local char- 
acteristics of the lattice site. The intramolecular force can be shown to satisfy a pressure 

Pk = Pkcl - ^ipl. (10) 


Equation 10 describes the equation of state (EoS) for the non-ideal EBE, and has both the ideal 
part and the non-ideal part. A real gas equation of state p = p{p) can be incorporated^^ if the 
effective mass is of the form 


■0fc 


' 2 {cjpk - Pk) 


( 11 ) 


In this case, phase segregation occurs consistent with the Maxwell construction when the tempera- 
ture is subcritical. The interaction potential is computed by approximating the spatial gradient on 
the lattice using a central finite difference scheme. 


While various equations of state can be incorporated into the EBE pseudopotential model through 
the use of Equation 11, the Carnahan-Starling EoS^^ has been widely used in the EBE due to its 
accuracy and numerical properties. Eor use in the EBE, the C-S EoS is often given in the pressure- 
density form as 


p = pRT 


1 + bp/4 + (6p/4)^ - {bp/ A)" 
(l-6p/4)3 


ap 


( 12 ) 


The C-S EoS is parameterized by the coefficients a, b, R, and the temperature T, but unlike 
some other EoS, has no free parameters. That is, a, b are completely specified by the critical point 
conditions. Since the EoS is valid in any unit system, for proper incorporation into the EBE with 
0 < p(xfc) < 1 the most common parameterization is to choose a = 1, 6 = 4, and R = 1. 

Accuracy of the pseudopotential model can be broadly characterized into two categories: (1) 
accuracy of the EoS in reproducing the commonly accepted or experimental data associated with a 
particular gas or liquid phase of a substance, and (2) accuracy of the EBE in converging the densities 
to those modeled by the EoS at a particular temperature and pressure. Both affect overall solution 
accuracy. 

A comparison with reference coexistence data for liquid oxygen was performed using the C-S 
EoS. The C-S EoS over-predicts the density of both the liquid and vapor phases and over-predicts 
the pressure. Since the simulation will be executed in isothermal and isochoric conditions, the 
C-S critical parameters were corrected to match the target parameters for the liquid at a specified 
temperature. The simulation results confirmed that the C-S EoS is only locally accurate, and still 
over-predicts the pressure slightly. This is a fundamental limitation in the use of a universal EoS. 
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3.2 Fluid- Wall Interaction 


Interaction with solid boundary conditions are specified using a supplementary wall wetting 
model that is similar to the pseudopotential force. First, the potential function near the wall must 
be modified in order to reconstruct the unknown gradient near the wall. Other than when explicitly 
specified, it is preferred that the fluid ignore the contribution of the wall. Local modification of the 
effective mass function is performed such that the gradient near the wall is identically zero and a 
supplemental wall interaction force is applied using a boundary condition function Sk equal to one 
if Xfc G d(f and zero otherwise. The parameter G^j is the wall interaction coefficient. 

The contact angle can be indirectly specified by Gw, where Gw < 0 implies a non- wetting fluid 
and Gw > 0 implies a wetting fluid. For Gw > 0, the interaction potential is attractive for the 
higher-density liquid phase, and thus the contact angle will approach zero. For Gw = 0, the wall 
behaves as the phase in which it is immersed and the contact angle will be small. For Gw < 0, the 
lower-density vapor phase will be attracted to the wall and the contact angle will be less than tt. In 
all cases, the wall interaction coefficient must be selected based upon the known wetting properties 
of the fluid under study. 


3.3 Dissipation and Stability 

In many implementations appearing in the literature, the stability is fundamentally limited by 
the kinematic viscosity of the modeled fluid, which is assumed to be uniform in the lattice. The 
development of instability as the density ratio is increased is usually the result of high velocities 
adjacent to the interface. Since the density of the vapor phase near the interface is very small, the 
high interface forces combined with the high-order anisotropy of the potential induce errors that 
lead to numerical instability.^ 

In the actual fluid, the invariant isothermal transport coefficient is the dynamic viscosity fi, and 
the kinematic viscosity is given by the relationship of the dynamic viscosity and the fluid density. 
Thus, the local relaxation rate can be modified to account for the variation in kinematic viscosity 
within the lattice. The local relaxation frequency is given by 


Wfc 


2 

6uk + 1 


(13) 


where is the local kinematic viscosity, Pk- The relaxation rates are specified on a per-cell 

basis using the method given in Equation 8. As the vapor phase density is much lower, the kinematic 
viscosity of the gas is much higher than that of the liquid. This model is consistent with the physics 
and substantially improves the stability near the interface. 

In certain applications where instabilities appear due to local density fluctuations or high Reynolds 
number flow, it may be appropriate to use a sub-grid turbulence scheme, such as the Smagorinsky 
model, to model eddy viscosity that exists at a scale smaller than can be resolved at a given lat- 
tice resolution.*^ The use of a turbulence model should be approached with caution as it must be 
tuned to accurately predict the physical behavior. However, the computational cost of the turbu- 
lence model with a coarse grid may be less than the fine grid required to resolve the detail, with 
both parameterizations yielding equivalent bulk fluid behavior. 
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4 IMPLEMENTATION AND EVALUATION 


In the present research, all of the LBE computation, data processing, and visualization is im- 
plemented directly in MATLAB. Extensive use of MATEAB’s multi-dimensional array operators 
allows many of the necessary operations, such as advection, collision, and the calculation of body 
forces, to be accomplished without the use of nested loops. By relying on MATEAB ’s internally 
optimized matrix libraries for much of the large-scale multiplication and division operations, a con- 
siderable advantage in computational efficiency and code simplicity is realized. All of the test cases 
detailed in this paper were run on Macintosh laptop computer with a 2.3 GHz Intel Core i7 CPU 
and 8 GB of 1.6 GHz DDR3 memory. 


4.1 Wall Wetting Effects 

Surface tension forces affecting fluid-wall interaction are important in the modeling of propellant 
management devices, such as traps and vanes. The behavior of the wall wetting model was inves- 
tigated using a special domain containing a variety of comer geometry. The domain size is 5 cm^ 
with a lattice dimension n = 128. The maximum relaxation frequency u = 1.998 giving a time 
step of f = 0.162 msec with no domain acceleration. A droplet of radius 20 lattice units was placed 
at the center of the domain and the sensitivity to was determined. 

The results are shown in Eigure 2. As the wall interaction potential becomes positive, the contact 
angle decreases and the fluid strongly adheres to acute corners. However, for G^j >0.1 and G^ < 
—0.5 a sharp increase in spurious velocities is noted, possibly due to the sharp gradients near the 
corners. In addition, for positive values of G^ the fluid is attracted to the wall, creating a nonuniform 
density near the wall due to the local compressibility error. 
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Figure 2. Sensitivity to wetting model parameter (left, = 0, right, Gm = —0.5) 


The use of Gw < 0, as expected, produces a hydrophobic wall and increases the contact angle 
to vr/2 or greater. No impacts to stability are observed for small values of Gw- Eor large negative 
values of Gw, the wall is strongly hydrophobic and a low-density layer will form between the droplet 
and the wall. 


4.2 Lateral Sloshing 


Lateral sloshing dynamics are a fundamental concern for control system stability in rocket-propelled 
space vehicles. Analytical solutions for lateral slosh motion in an axisymmetric container can be 
derived for a variety of conditions; the simplest case describes a container in a high-g condition 
where the free surface height h is greater than the diameter, h > 2a where a is the radius.^ In this 
case, the bottom geometry can be ignored, the contact angle is near tt/2, and the first mode is a 
planar axisymmetric motion of the free surface along a fixed node line. In most cases, this simple 
linear mechanical analog of the fluid motion can be used to approximate the forces and moments on 
the tank. Importantly for many analyses, the mechanical model usually predicts the frequency to a 
reasonable level of accuracy. 

As the Bond number decreases below about Bo = 1000, the contact angle assumption must be 
modified and extensions to the basic potential theory must be employed. In addition, the presence 
of hemispherical domes or other special geometry necessitates corrections to the parameters of the 
equivalent mechanical model when the fluid level is near the upper or lower dome.^ However, 
analytical models can still reasonably predict the sloshing parameters for Bond numbers as low as 
Bo = 10.^^ 

The lateral sloshing test case is constructed in a closed domain taken as the cross section of an 
axisymmetric cylindrical tank having hemispherical ends. The tank radius is a = 0.15 m and the 
barrel section length is 0.2 m. The tank is under a static acceleration field of g = O.OOl^fQ and the 
simulated fluid lamina has a total mass of 0.1061 kg. The lattice size is n = 346 with a time step 
of 8t = 2.2 msec. Assuming that the surface tension is approximately a = 0.0122 N/m, the Bond 
number Bo ss 20. At this condition, the Bond number is low enough to be simulated by the LBE 
and near the lower limit of the analytic prediction of the first mode natural frequency. Since both 
analytic and numerical models can be used to predict the sloshing response, this presents a useful 
verification opportunity. 

The first lateral sloshing mode is excited by specifying an initial condition with the acceleration 
rotated +15° from the vertical orientation. The system is then allowed to reach a quasi-steady state 
condition with a total runtime of f = 240 seconds. To simulate free decay, the acceleration is 
returned to an angle aligned with the tank symmetry axis. The free decay is simulated for f = 60 
seconds to capture several cycles of the very long-period fundamental lateral slosh wave. The 
density and lattice velocity field after t = 9.8 seconds are shown in Figure 3. The time history of the 
force decay and spectral analysis of the domain response are shown in Figure 4. The fundamental 
mode is clearly visible (and appears in the frequency response) but is not dominant in frequency 
spectrum compared with a higher-frequency traveling-wave type free-surface mode. It is notable 
that the axial component of the response has a noticeable oscillation at one half the higher lateral 
mode frequency. These higher modes are probably an artifact of the initial conditions, which have 
not completely decayed, and the high-frequency response of the liquid to an instantaneous shift in 
the acceleration vector. 

At a Bond number near Bo = 20, the high-g analytic solution for a right circular cylinder is a 
reasonable prediction for the fundamental antisymmetric mode frequency; 

tanh f h}il\ (^ 14 ) 

a \ a J 

where are the roots of the eigenvalue equation dJ\{\r / a) / dr = 0 for r = a and J\ is a Bessel 
function of the first kind. At low Bond number and for a liquid depth of h > 3a, the empirical 
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Figure 3. Lateral sloshing response after 9.8 seconds 


relation* for the first mode 

w? = 2.592-^(l + 0.798So) (15) 

pa-^ 

correlates well with experiment; the frequency tends to increase slightly as the Bond number de- 
creases. 

The frequency also has an exact solution if the contact angle 6c = -n / 2? In this case, 


= 


2 \2 
n 


a 


pa-^ 




An tanh 


Xrih 


(16) 


Finally, the analytical solution corrected for the bottom geometry'^ can give an estimate of the 
natural frequency at this condition. The predicted and simulated frequencies are summarized in 
Table 1. The LBE result agrees very well with the analytical predictions. 


Method 

Mode n = 1 frequency, Hz 

Equation 14 

0.0552 

Equation 15 

0.0603 

Equation 16 

0.0595 

Dome corrected (reference'^) 

0.0551 

EBE (present result) 

0.0552 


Table 1. Predicted sloshing natural frequencies, ~50% fill 


* A typographical error appears in Equation 15 in the original reference.^ 


10 





Figure 4. Time history of tank forces (left) and frequency response (right) 


5 DISCUSSION AND CONCLUSIONS 

It has been shown that the lattice Boltzmann method has the potential to support modeling of low 
Bond number flows with applications to cryogenic spacecraft propellant dynamics. Some compu- 
tational and analytical advantages of the LBE have been highlighted, including data locality and 
explicit incorporation of an isothermal equation of state. 

However, the lattice Boltzmann approach requires further study to increase its level of maturity. 
The LBE in the present implementation has some drawbacks that limit its scope in application to 
the target flow regimes. These include numerical instability for steady flow velocities exceeding 
about 10% of the lattice speed of sound, implicit dependence of the surface tension on the lattice 
parameters, and a limited ability to incorporate convective effects. 

Many of these shortcomings have been addressed in the literature and can be resolved through ad- 
ditional development. Eirst, the pseudopotential method provides an ability to incorporate arbitrary 
EoS into the model. ^ ' Eor certain species such as oxygen, other parameterized EoS may yield better 
accuracy over a wide range of temperatures.'^ The use of adaptive time stepping is straightforward 
in the EBE and has shown some promise;'^ while the Mach number of the flow is a function of 
the adaptive parameter, this allows larger steps to be taken when the flow is quiescent and u <C 1. 
In this case, the accuracy may even be improved because the truncation error in the higher-order 

velocity terms can be kept fixed. Einally, various multi-phase thermal lattice Boltzmann techniques 
19 20 

are emerging. ’ 

In summary, the EBE shows promise in application to microgravity fluid flow regimes. The im- 
provement of accuracy and stability using the aforementioned methods and the continued validation 
of the results using theory and test data are the subjects of future study. 
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